%% EXERCISE 31: Cantilever beam - students's version

clear all
close all
clc

format short
addpath('ex31_toolbox_local')


%% Partition of unity option

%% LME options (Machine precision)
optLME.dim   = 2;
optLME.verb  = 0; %0:off
optLME.grad  = 1;           % Computation of the Gradient 0:OFF 1:ON
optLME.hess  = 0;           % Computation of the Hessian  0:OFF 1:ON
optLME.TolNR = 1.e-10;
optLME.knn   = 0;
optLME.Tol0  = 1.e-06;
optLME.gamma = 4.8;      % value of gamma to compute LME

% Nodal Integration options
%cubature order: 1 (1 GPts), 2 (3 GPts), 4 (6 GPts), 5 (7 GPts), 6 (12 GPts)
optGL.orderGL = 4;

%% Matrial Parameters
E  = 1000;     % Young modulus
nu = 0.3;     % Poisson coefficient
parameters.E  = E;
parameters.nu = nu;

%% 
Lx = 3;
Ly = 3;

% Node Points
Ny = 4;
Nx = 4;

h  = Lx /(Nx-1);
optLME.spacing = h;

% This line is important
center  = 0.5*[Lx Ly];
x_nodes = UniformGrid2D(Nx, Ny, Lx, Ly, center);

nPts  = length(x_nodes);

ids   = (1:nPts);
id_bd = [ids(x_nodes(:,1)==0) ids(x_nodes(:,2)==0) ...
    ids(x_nodes(:,1)==Lx) ids(x_nodes(:,2)==Ly)];
    id_bd = unique(id_bd);

% Sample points and gauss weights
[x_samples w_samples] = MakeGLSamples2D(x_nodes, optGL);

sPts = length(x_samples);


%   Patch transform
%A = [1          0;...
%     0          1         ]; % No scaling
 A = [1          sqrt(2)/2;...
      -sqrt(2)/2  1/2         ]; % No scaling

parameters.ind_Dirichlet = id_bd;
parameters.linear_transform = A;
%% Nodes thermalization, nodes and samples adjacency structures are computed
[beta_n range_n] = NodalThermalization(x_nodes, optLME);
optLME.beta   = optLME.gamma/(h*h);
optLME.beta_n = beta_n;
optLME.range_n= range_n;
  
%% The basis functions are computed
% adjacency structure with the nearest neighbors nodes to each sample point
% nodal shape parameter
s_near = SamplesAdjacency(x_nodes,x_samples,range_n);
  
% Local-max entropy basis functions computation
optLME.s_near = s_near;
outLME = wrapper_lme(x_nodes,x_samples,optLME);
  
p_samp  = outLME.p_samp;
dp_samp = outLME.dp_samp;
fprintf(1,'cputime LME   : %4.3f\n', cputime);

%% The displacement field is computed
[u_h u_exact] = ex31_SolveSystem(dp_samp,s_near,x_nodes,x_samples,w_samples,parameters,optLME);
%%
figure(1);
cla;
scatter(x_nodes(:,1),x_nodes(:,2),'g'); hold on;
quiver(x_nodes(:,1),x_nodes(:,2),u_exact(1:2:end-1),u_exact(2:2:end),'r');
title('Exact solution');
figure(2);
cla;
scatter(x_nodes(:,1),x_nodes(:,2),'g'); hold on;
quiver(x_nodes(:,1),x_nodes(:,2),u_h(1:2:end-1),u_h(2:2:end),'b');
title('Aproximated solution');

fprintf(1,'Infinite norm of the error: %d \n',max(u_h-u_exact));
